import numpy as np
from EOS_collection import NRTL_MATRIX, antoine, gamma_cal, G_ex_MATRIX
import Mixing_rules as MR
from scipy.optimize import minimize

# 常数
R = 8.314
atm = 101.325
Pa = atm
alpha = 0.3

# 基础参数
a = np.array([1.201, -0.317])
b = np.array([1.001, 0.995])
T = 365.54
z = np.array([[0.1, 0.9]])
z = z.reshape(-1, 1)
x = z
y = z
LtoF = 0.5

A = np.array([7.19736, 7.07404])
B = np.array([1574.99, 1657.46])
C = np.array([-34.29, -46.13])
Pvs = antoine(A, B, C, T)

# 计算活度系数
tau_unknown = a + b / T
tau = np.zeros((2, 2))
tau[np.where(np.eye(2)==0)] = tau_unknown.flatten() #除了主对角线，其余依次赋值

def Fraction(LtoF, x, y, z, P, T):
    # 计算液相逸度系数
    sumxy = 1
    LtoF_new = LtoF

    while abs(sumxy) > 1e-8:
        
        LtoF = LtoF_new

        gamma = NRTL_MATRIX(x.T.reshape(-1, 1), tau, alpha).flatten() #将x变成列向量

        k = Pvs / Pa * gamma
        x = z / (k.reshape(-1, 1) + LtoF * (1 - k.reshape(-1, 1)))
        y = k.reshape(-1, 1) * x

        sumxy = np.sum((k - 1) * z.flatten() / (k + LtoF * (1 - k)))
        term1 = sumxy
        term2 = np.sum(z.flatten() * (k - 1) ** 2 / (k + (1 - k) * LtoF) ** 2)
        LtoF_new = LtoF - term1 / term2
        # print(f"term1 = {sumxy}")
        # print(f"term2 = {term2}")

    print(f"液相组成 = {x}")
    print(f"气相组成 = {y}")
    return LtoF_new

LtoF = Fraction(LtoF, x, y, z, Pa, T)
print(f"液相分率是{LtoF}")